#include "../include/Distributions.hpp"
#include "../include/Density.hpp"
#include "../include/Kernel.hpp"
#include "../include/Particle.hpp"
#include "../include/Resample.hpp"
#include "../include/utils.hpp"
#include "../include/Models/Likelihood.hpp"
#include "../include/Misc/EP.hpp"
#include "../include/Misc/Laplace.hpp"
#include "../include/Models/ProbitMS.hpp"
#include "../include/Models/ProbitMS2.hpp"
#define NU 200

extern "C"{

void IBIS(int *M, int *n, int *p, double *Y, double *pX,double *res, double *W)
{
	mat X=vect2mat<double>(pX,n[0],p[0]);
//	cout << X;
	//Density::Density_Probit *DP= new Density::Density_Probit(Y,X);
	mat m(p[0],1);
	mat foo(p[0],p[0]);
	m.zeros();
	foo.fill(0.01);
	mat s=diagmat(foo);
	
	mat bar(p[0],p[0]);
	bar.fill(1);
	mat S=diagmat(bar);

//	EP fEP;
//	fEP(X,Y,0.1,3);

	Laplace fL;
	fL(X,Y,m,S);	

	
	mat Sig=fL.Get_Sig();	
	mat Mu=fL.Get_Mu();

	//	cout << s;
	Distribution::Uniform *U2;
	Distribution::LogitProp *LP1; 
//	
	
	Probit_MS *PMS;
	Kernel_MHID<Density::GeomBridge,Distribution::LogitProp> *K;
	//Kernel_MH<Density::GeomBridge,Distribution::Student> *K= new Kernel_MH<Density::GeomBridge,Distribution::Student>(DP,G)
	//Resample::SysResample *R=new Resample::SysResample;
	Resample::Multinomial *R;
	//Particle<Kernel_MH<Density::IBIS,Distribution::Gaussian>,Resample::Multinomial,Density::IBIS> P(K,R,DP,M[0]);
	//Particle<Kernel_MH<Density::Tempering,Distribution::Gaussian>,Resample::SysResample,Density::Tempering> P(K,R,DP,M[0]);
//	Particle<Kernel_MH<Density::GeomBridge,Distribution::Student>,Resample::Multinomial,Density::GeomBridge> P(K,R,DP,M[0],F);
//	cout << "hello";

	//mat gamma=U2->r(1);
	//cout << gamma ;
//for( int i=0;i<1000;i++){double t=PMS->Likelihood(gamma);}

	Particle<Kernel_MHID<Density::GeomBridge,Distribution::LogitProp>,Resample::Multinomial,Density::GeomBridge> *P;
	//http://openmp.org/wp/resources
	int numCPU = sysconf( _SC_NPROCESSORS_ONLN );
	//numCPU=1;
	cout << "cpu  " << numCPU;
	mat A((int)(M[0]/numCPU),p[0]-1);
	mat w((int)(M[0]/numCPU),1);
	vector<mat> weights(numCPU,w);
	vector<mat> Pi(numCPU,A);

	#pragma omp parallel for firstprivate(R,K,U2,LP1,PMS,P), shared(Pi,weights,X,Y)
	for(int i=0;i<numCPU;i++)
	{
		U2=new Distribution::Uniform(p[0]-1,0,1);
		LP1=new Distribution::LogitProp(p[0]-1);
		PMS=new Probit_MS(Y,X,U2,S);
		K= new Kernel_MHID<Density::GeomBridge,Distribution::LogitProp>(PMS,LP1);
		R=new Resample::Multinomial;
		P=new Particle<Kernel_MHID<Density::GeomBridge,Distribution::LogitProp>,Resample::Multinomial,Density::GeomBridge>(K,R,PMS,M[0]/numCPU,U2);
		P->Filter();
		Pi[i]=P->Get_theta();
		weights[i]=P->Get_W();
	}

	mat Wi=vectortomat(weights,numCPU);
	mat theta=vectortomat(Pi,numCPU);
	
	
//	cout << DP->Density(A.row(0).t(),0,n[0]-1);	
//	cout << A;
//	cout << "hello";
//cout << M[0] << p[0];	
	mat2vect<double>(theta,M[0],p[0]-1,res);
	mat2vect<double>(Wi,M[0],1,W);
/*
	delete LP1;
//	delete DP;
	delete U2;
	delete R;
	delete K;
*/	
}

}
